# Dickstein, Ho, and Mark (2023)
# This script is called by EA02*
# It estimates small group surplus for the bootstraped results.

# * # * # * # * # * # * # * # * #
# DERIVING V ESTIMATES          #
# * # * # * # * # * # * # * # * #

# Subsetting exploded data to the relevant sample
ExplodedData_SG <-  mod_exploded_data[best_guess_ra %in% ravec]
ExplodedData_SG <- ExplodedData_SG[year %in% yrvec]

# Setting SG market demand parameter values
SGNamedBetaList <- list(
  beta1 = as.numeric(SGMarketEstimates[1, 1:SGMarSE[1]]),
  beta2 = as.numeric(SGMarketEstimates[1, (1 + SGMarSE[1]):SGMarSE[2]]),
  beta3 = as.numeric(SGMarketEstimates[1, (1 + SGMarSE[2]):SGMarSE[3]]),
  beta0 = as.numeric(SGMarketEstimates[1, (1 + SGMarSE[3]):SGMarSE[4]])
)

# Setting SG market demand parameter names
names(SGNamedBetaList[[1]]) <- names(SGMarketEstimates)[1:SGMarSE[1]]
names(SGNamedBetaList[[2]]) <- names(SGMarketEstimates)[(1 + SGMarSE[1]):SGMarSE[2]]
names(SGNamedBetaList[[3]]) <- names(SGMarketEstimates)[(1 + SGMarSE[2]):SGMarSE[3]]
names(SGNamedBetaList[[4]]) <- names(SGMarketEstimates)[(1 + SGMarSE[3]):SGMarSE[4]]

# Construct alpha, omega, psi, fixed effect term in the Exploded Dataset:
ExplodedData_SG$alpha <- exp(as.matrix(ExplodedData_SG[, names(SGNamedBetaList$beta1), with = F]) %*% SGNamedBetaList$beta1)
ExplodedData_SG$omega <- exp(as.matrix(ExplodedData_SG[, names(SGNamedBetaList$beta2), with = F]) %*% SGNamedBetaList$beta2)
ExplodedData_SG$psi <- exp(as.matrix(ExplodedData_SG[, names(SGNamedBetaList$beta3), with = F]) %*% SGNamedBetaList$beta3)
ExplodedData_SG$fixedeffect_beta0 <- as.matrix(ExplodedData_SG[, names(SGNamedBetaList$beta0), with = F]) %*% SGNamedBetaList$beta0

# Generating Household Level Price: 
## no tiering
ExplodedData_SG[, price_hh_notiering := (grossprem / 100) * (1 - nu) * (1 - (best_guess_sumrate / 100)), ]
## tiering
ExplodedData_SG[, totalprem_contract := sum(grossprem), by = c("constructed_plan_year", "contractnum", "year")]
ExplodedData_SG[, totalfamtypescore_contract := sum(famtype_score), by = c("constructed_plan_year", "contractnum", "year")]
ExplodedData_SG[, grossprem_tiered := famtype_score * totalprem_contract / totalfamtypescore_contract, ]
ExplodedData_SG[, price_hh := (grossprem_tiered / 100) * (1 - nu) * (1 - (best_guess_sumrate / 100)), ]
### testing for accuracy
print(paste0(
  "The sum of no tiered prices is ", 
  sum(ExplodedData_SG$price_hh_notiering) / 1000000, 
  " million. The sum of tiered prices is ", 
  sum(ExplodedData_SG$price_hh) / 1000000,
  " million. They should be close but not equal."
)) 

# Generating V: 
ExplodedData_SG[, V_notiering := x_av + .5 * omega * (x_av ^ 2) - (alpha - psi) * price_hh_notiering + fixedeffect_beta0, ]
ExplodedData_SG[, V := x_av + .5 * omega * (x_av ^ 2) - (alpha - psi) * price_hh + fixedeffect_beta0, ]

# * # * # * # * # * # * # * # * #
# SAVING DATA FOR PROCESS       #
# * # * # * # * # * # * # * # * #
ExplodedData_SG <- merge(ExplodedData_SG, 
  SubDat, 
  by.x = c("subscriberid", "year"), 
  by.y = c("subscriberid", "year"), 
  all.x = T, 
  suffixes = c("", ".subdat"))

ExplodedData_SG[, choice := as.numeric(chosen_plan_id == constructed_plan_year), ]
nchoice <- ExplodedData_SG[, (nchoice = .N), by = c("subscriberid", "year")]
print(paste0(
  "There are ", 
  sum(nchoice$nchoice != 1), 
  " household with number of choices not equal to 1"
))

# Saving
ExplodedData_SG[, bronze := as.numeric(x_av == .6), ]
ExplodedData_SG[, silver := as.numeric(x_av == .7 | x_av == .73 | x_av == .87), ]
ExplodedData_SG[, gold := as.numeric(x_av == .8), ]
ExplodedData_SG[, platinum := as.numeric(x_av == .9), ]

varstosave <- c("subscriberid", "year", "constructed_plan_year", "choice", 
  "x_av", "group_size", "sum_concurrent_risk", "max_concurrent_risk", 
  "age", "family_size", "familytype", "nspouse", "ndeps", 
  "bronze", "silver", "gold", "platinum", "price_hh_notiering", "grossprem", "grossprem_tiered", "price_hh",
  "alpha", "psi", "omega", "fixedeffect_beta0", "V", "V_notiering")
sg_surplus <- ExplodedData_SG[, varstosave, with = F]

sg_surplus <- sg_surplus[, hhid := .GRP, by = c("subscriberid", "year")]
